

namespace Tuna {

template<class Tprec, int Dim>
inline bool Simplec<Tprec, Dim>::calcCoefficients2D()
{
  Tprec ce, cw;
  Tprec cn, cs;  
  aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;
    

  int pausa;

  for (int i =  bi; i <= ei; ++i) 
    for (int j = bj; j <= ej; ++j) 
      {
	ce = u(i  ,j) * dy;
	cw = u(i-1,j) * dy;
	cn = v(i,j  ) * dx;
	cs = v(i,j-1) * dx;
	
	aE (i,j) = du(i,j) * dy;
	aW (i,j) = du(i-1,j) * dy;
	aN (i,j) = dv(i,j) * dx;
	aS (i,j) = dv(i,j-1) * dx;
	aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j);
	sp(i,j) = cw - ce + cs - cn;
      }

  return 0;
}

template<class Tprec, int Dim>
inline bool Simplec<Tprec, Dim>::calcCoefficients3D()
{
    Tprec dyz = dy * dz;
    Tprec dxz = dx * dz;
    Tprec dxy = dx * dy;
    Tprec ce, cw;
    Tprec cn, cs;
    Tprec cf, cb;

    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0;
    aP = 0.0; sp = 0.0;
    
    for (int i = bi; i <= ei; ++i) 
	for (int j = bj; j <= ej; ++j) 	
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = u(i  ,j, k) * dyz;
		cw = u(i-1,j, k) * dyz;
		cn = v(i,j,k  ) * dxz;
		cs = v(i,j-1,k) * dxz;
		cf = w(i,j,k  ) * dxy;
		cb = w(i,j,k-1) * dxy;
				
		aE (i,j,k) = du(i,j,k) * dyz;
		aW (i,j,k) = du(i-1,j,k) * dyz;
		aN (i,j,k) = dv(i,j,k) * dxz;
		aS (i,j,k) = dv(i,j-1,k) * dxz;
		aF (i,j,k) = dw(i,j,k) * dxy;
		aB (i,j,k) = dw(i,j,k-1) * dxy;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + aS (i,j,k)
		    + aF (i,j,k) + aB (i,j,k);
		sp(i,j,k) = cw - ce + cs - cn + cb - cf;
	    }

    return 0;
}

template<class Tprec, int Dim>
bool Simplec<Tprec, Dim>::correction() {

///// ----- Velocity correction -----------------------
// 
// OJO CON LAS COTAS, SON DIFERENTES PARA u, v, w y p!!!!
//
    Range Iu( u.lbound(firstDim) + 1, u.ubound(firstDim) - 1  );
    Range Ju( u.lbound(secondDim) + 1, u.ubound(secondDim) - 1 );
    Range Iv( v.lbound(firstDim) + 1, v.ubound(firstDim) - 1  );
    Range Jv( v.lbound(secondDim) + 1, v.ubound(secondDim) - 1 );
    
    if (Dim == 2) {
	u(Iu, Ju) = u(Iu, Ju) + 
	    du(Iu, Ju) * ( phi(Iu, Ju) - phi(Iu+1, Ju) );
	
	v(Iv, Jv) = v(Iv, Jv) + 
	    dv(Iv, Jv) * ( phi(Iv, Jv) - phi(Iv, Jv+1) );
    }
    if (Dim == 3) {
	Range Iw( w.lbound(firstDim) + 1, w.ubound(firstDim) - 1  );
	Range Jw( w.lbound(secondDim) + 1, w.ubound(secondDim) - 1 );
	Range Kw( w.lbound(thirdDim) + 1, w.ubound(thirdDim) - 1 );
	Range Ku( u.lbound(thirdDim) + 1, u.ubound(thirdDim) - 1 );
	Range Kv( v.lbound(thirdDim) + 1, v.ubound(thirdDim) - 1 );
	
	u(Iu, Ju, Ku) = u(Iu, Ju, Ku) + 
	    du(Iu, Ju, Ku) * ( phi(Iu, Ju, Ku) - phi(Iu+1, Ju, Ku) );
	
	v(Iv, Jv, Kv) = v(Iv, Jv, Kv) + 
	    dv(Iv, Jv, Kv) * ( phi(Iv, Jv, Kv) - phi(Iv, Jv+1, Kv) );
	
	w(Iw, Jw, Kw) = w(Iw, Jw, Kw) + 
	    dw(Iw, Jw, Kw) * ( phi(Iw, Jw, Kw) - phi(Iw, Jw, Kw+1) );
    }

///// ----- Pressure correction -----------------------
    intTinyArray_t lower_bounds, upper_bounds;
    if (Dim == 2) { 
	lower_bounds = phi.lbound(firstDim), phi.lbound(secondDim);
	upper_bounds = phi.ubound(firstDim), phi.ubound(secondDim);
    }
    if (Dim == 3) {
	lower_bounds = phi.lbound(firstDim), 
	    phi.lbound(secondDim),
	    phi.lbound(thirdDim);
	upper_bounds = phi.ubound(firstDim),
	    phi.ubound(secondDim),
	    phi.ubound(thirdDim); 
    }
    RectDomain<Dim> domain(lower_bounds, upper_bounds);
//
// phi <-- variable of equation, in this case pressure correction p'
// p   <-- pressure
//
    Tprec alpha = 1.000; // this will be used for under-relaxation...


    p(domain) = p(domain) + alpha * phi(domain);




    phi = 0;



    return 0;
}

} // Tuna namespace
